---
title: "Subgroup analysis: risk talk"
author: "[redacted]"
date: "2025-04-07"
output: html_document
---

```{r}
```

#####LOAD PACKAGES AND DATA

```{r}
#Load packages
library(tidyverse)
library(cjoint)
library(cregg)

conjoint_rp <- cjoint::read.qualtrics(
  "data_clean.csv", 
  covariates = c(
                  "risk_prob_1_fipl_1", 
                  "risk_prob_1_fipl_2", 
                  "risk_prob_1_fipl_3", 
                  "risk_prob_1_fipl_4", 
                  "risk_sev_1_fipl_1", 
                  "risk_sev_1_fipl_2", 
                  "risk_sev_1_fipl_3", 
                  "risk_sev_1_fipl_4",
                  "risk_prob_2_filp_1", 
                  "risk_prob_2_filp_2", 
                  "risk_prob_2_filp_3", 
                  "risk_prob_2_filp_4",
                  "risk_sev_2_filp_1", 
                  "risk_sev_2_filp_2", 
                  "risk_sev_2_filp_3", 
                  "risk_sev_2_filp_4",
                  "risk_prob_3_fpil_1", 
                  "risk_prob_3_fpil_2", 
                  "risk_prob_3_fpil_3", 
                  "risk_prob_3_fpil_4",
                  "risk_sev_3_fpil_1", 
                  "risk_sev_3_fpil_2", 
                  "risk_sev_3_fpil_3", 
                  "risk_sev_3_fpil_4",
                  "risk_prob_4_fpli_1", 
                  "risk_prob_4_fpli_2", 
                  "risk_prob_4_fpli_3", 
                  "risk_prob_4_fpli_4",
                  "risk_sev_4_fpli_1", 
                  "risk_sev_4_fpli_2", 
                  "risk_sev_4_fpli_3", 
                  "risk_sev_4_fpli_4",
                  "risk_prob_5_flip_1",         
                  "risk_prob_5_flip_2", 
                  "risk_prob_5_flip_3", 
                  "risk_prob_5_flip_4",
                  "risk_sev_5_flip_1", 
                  "risk_sev_5_flip_2", 
                  "risk_sev_5_flip_3", 
                  "risk_sev_5_flip_4",
                  "risk_prob_6_flpi_1",
                  "risk_prob_6_flpi_2", 
                  "risk_prob_6_flpi_3", 
                  "risk_prob_6_flpi_4",
                  "risk_sev_6_flpi_1", 
                  "risk_sev_6_flpi_2", 
                  "risk_sev_6_flpi_3", 
                  "risk_sev_6_flpi_4",
                  "risk_prob_7_ifpl_1",
                  "risk_prob_7_ifpl_2", 
                  "risk_prob_7_ifpl_3", 
                  "risk_prob_7_ifpl_4",
                  "risk_sev_7_ifpl_1", 
                  "risk_sev_7_ifpl_2", 
                  "risk_sev_7_ifpl_3", 
                  "risk_sev_7_ifpl_4",
                  "risk_prob_8_iflp_1",
                  "risk_prob_8_iflp_2", 
                  "risk_prob_8_iflp_3", 
                  "risk_prob_8_iflp_4",
                  "risk_sev_8_iflp_1", 
                  "risk_sev_8_iflp_2", 
                  "risk_sev_8_iflp_3", 
                  "risk_sev_8_iflp_4",
                  "risk_prob_9_ipfl_1",
                  "risk_prob_9_ipfl_2", 
                  "risk_prob_9_ipfl_3", 
                  "risk_prob_9_ipfl_4",
                  "risk_sev_9_ipfl_1", 
                  "risk_sev_9_ipfl_2", 
                  "risk_sev_9_ipfl_3", 
                  "risk_sev_9_ipfl_4",
                  "risk_prob_10_iplf_1",
                  "risk_prob_10_iplf_2", 
                  "risk_prob_10_iplf_3", 
                  "risk_prob_10_iplf_4",
                  "risk_sev_10_iplf_1", 
                  "risk_sev_10_iplf_2", 
                  "risk_sev_10_iplf_3", 
                  "risk_sev_10_iplf_4",
                  "risk_prob_11_ilfp_1",
                  "risk_prob_11_ilfp_2", 
                  "risk_prob_11_ilfp_3", 
                  "risk_prob_11_ilfp_4",
                  "risk_sev_11_ilfp_1", 
                  "risk_sev_11_ilfp_2", 
                  "risk_sev_11_ilfp_3", 
                  "risk_sev_11_ilfp_4",
                  "risk_prob_12_ilpf_1",
                  "risk_prob_12_ilpf_2", 
                  "risk_prob_12_ilpf_3", 
                  "risk_prob_12_ilpf_4",
                  "risk_sev_12_ilpf_1", 
                  "risk_sev_12_ilpf_2", 
                  "risk_sev_12_ilpf_3", 
                  "risk_sev_12_ilpf_4",
                  "risk_prob_13_pfil_1",
                  "risk_prob_13_pfil_2", 
                  "risk_prob_13_pfil_3", 
                  "risk_prob_13_pfil_4",
                  "risk_sev_13_pfil_1", 
                  "risk_sev_13_pfil_2", 
                  "risk_sev_13_pfil_3", 
                  "risk_sev_13_pfil_4",
                  "risk_prob_14_pfli_1",
                  "risk_prob_14_pfli_2", 
                  "risk_prob_14_pfli_3", 
                  "risk_prob_14_pfli_4",
                  "risk_sev_14_pfli_1", 
                  "risk_sev_14_pfli_2", 
                  "risk_sev_14_pfli_3", 
                  "risk_sev_14_pfli_4",
                  "risk_prob_15_pifl_1",
                  "risk_prob_15_pifl_2", 
                  "risk_prob_15_pifl_3", 
                  "risk_prob_15_pifl_4",
                  "risk_sev_15_pifl_1", 
                  "risk_sev_15_pifl_2", 
                  "risk_sev_15_pifl_3", 
                  "risk_sev_15_pifl_4",
                  "risk_prob_16_pilf_1",
                  "risk_prob_16_pilf_2", 
                  "risk_prob_16_pilf_3", 
                  "risk_prob_16_pilf_4",
                  "risk_sev_16_pilf_1", 
                  "risk_sev_16_pilf_2", 
                  "risk_sev_16_pilf_3", 
                  "risk_sev_16_pilf_4",
                  "risk_prob_17_plfi_1",
                  "risk_prob_17_plfi_2", 
                  "risk_prob_17_plfi_3", 
                  "risk_prob_17_plfi_4",
                  "risk_sev_17_plfi_1", 
                  "risk_sev_17_plfi_2", 
                  "risk_sev_17_plfi_3", 
                  "risk_sev_17_plfi_4",
                  "risk_prob_18_plif_1",
                  "risk_prob_18_plif_2", 
                  "risk_prob_18_plif_3", 
                  "risk_prob_18_plif_4",
                  "risk_sev_18_plif_1", 
                  "risk_sev_18_plif_2", 
                  "risk_sev_18_plif_3", 
                  "risk_sev_18_plif_4",
                  "risk_prob_19_lfip_1",
                  "risk_prob_19_lfip_2", 
                  "risk_prob_19_lfip_3", 
                  "risk_prob_19_lfip_4",
                  "risk_sev_19_lfip_1", 
                  "risk_sev_19_lfip_2", 
                  "risk_sev_19_lfip_3", 
                  "risk_sev_19_lfip_4",
                  "risk_prob_20_lfpi_1",
                  "risk_prob_20_lfpi_2", 
                  "risk_prob_20_lfpi_3", 
                  "risk_prob_20_lfpi_4",
                  "risk_sev_20_lfpi_1", 
                  "risk_sev_20_lfpi_2", 
                  "risk_sev_20_lfpi_3", 
                  "risk_sev_20_lfpi_4",
                  "risk_prob_21_lifp_1",
                  "risk_prob_21_lifp_2", 
                  "risk_prob_21_lifp_3", 
                  "risk_prob_21_lifp_4",
                  "risk_sev_21_lifp_1", 
                  "risk_sev_21_lifp_2", 
                  "risk_sev_21_lifp_3", 
                  "risk_sev_21_lifp_4",
                  "risk_prob_22_lipf_1",
                  "risk_prob_22_lipf_2", 
                  "risk_prob_22_lipf_3", 
                  "risk_prob_22_lipf_4",
                  "risk_sev_22_lipf_1", 
                  "risk_sev_22_lipf_2", 
                  "risk_sev_22_lipf_3", 
                  "risk_sev_22_lipf_4",
                  "risk_prob_23_lpfi_1",
                  "risk_prob_23_lpfi_2", 
                  "risk_prob_23_lpfi_3", 
                  "risk_prob_23_lpfi_4",
                  "risk_sev_23_lpfi_1", 
                  "risk_sev_23_lpfi_2", 
                  "risk_sev_23_lpfi_3", 
                  "risk_sev_23_lpfi_4",
                  "risk_prob_24_lpif_1",
                  "risk_prob_24_lpif_2", 
                  "risk_prob_24_lpif_3", 
                  "risk_prob_24_lpif_4",
                  "risk_sev_24_lpif_1", 
                  "risk_sev_24_lpif_2", 
                  "risk_sev_24_lpif_3", 
                  "risk_sev_24_lpif_4"),  
  respondentID = "ResponseId",
  letter = "F",
  new.format = TRUE,
  responses = paste0("conjoint_", c("intro", 1:8))
)

#For some reason, R keeps changing the name of the variable "risk_prob_1_fipl_2" to "respondent". I present the ingenious solution here:
colnames(conjoint_rp)[colnames(conjoint_rp) == "respondent"] <- "risk_prob_1_fipl_2"

names(conjoint_rp)


```

#RECODE COVARIATE

```{r}

###RISK PROBABILITY
# Create empty variables for the original questions in the conjoint dataset
conjoint_rp$rprob_flooding <- NA
conjoint_rp$rprob_infectious <- NA
conjoint_rp$rprob_pesticides <- NA
conjoint_rp$rprob_lifestyle <- NA
conjoint_rp$rs_flooding <- NA
conjoint_rp$rs_infectious <- NA
conjoint_rp$rs_pesticides <- NA
conjoint_rp$rs_lifestyle <- NA

# Map each block to its order
question_order <- c(
  'fipl', 'filp', 'fpil', 'fpli', 'flip', 'flpi',
  'ifpl', 'iflp', 'ipfl', 'iplf', 'ilfp', 'ilpf',
  'pfil', 'pfli', 'pifl', 'pilf', 'plfi', 'plif',
  'lfip', 'lfpi', 'lifp', 'lipf', 'lpfi', 'lpif'
)

# Convert all risk probability and severity columns to numeric
for (block in 1:24) {
  for (position in 1:4) {
    risk_prob_var <- paste0("risk_prob_", block, "_", question_order[block], "_", position)
    risk_sev_var <- paste0("risk_sev_", block, "_", question_order[block], "_", position)
    
    # Check if the columns exist and have data
    if (risk_prob_var %in% colnames(conjoint_rp)) {
      if (length(conjoint_rp[[risk_prob_var]]) == 16020) {  # Check if it has the expected length
        conjoint_rp[[risk_prob_var]] <- as.numeric(as.character(conjoint_rp[[risk_prob_var]]))
      } else {
        print(paste("Warning: Column", risk_prob_var, "does not have 16020 rows"))
      }
    } else {
      print(paste("Warning: Column", risk_prob_var, "does not exist in the dataset"))
    }

    if (risk_sev_var %in% colnames(conjoint_rp)) {
      if (length(conjoint_rp[[risk_sev_var]]) == 16020) {
        conjoint_rp[[risk_sev_var]] <- as.numeric(as.character(conjoint_rp[[risk_sev_var]]))
      } else {
        print(paste("Warning: Column", risk_sev_var, "does not have 16020 rows"))
      }
    } else {
      print(paste("Warning: Column", risk_sev_var, "does not exist in the dataset"))
    }
  }
}

colnames(conjoint_rp)


# Convert all risk probability and severity columns to numeric
for (block in 1:24) {
  for (position in 1:4) {
    risk_prob_var <- paste0("risk_prob_", block, "_", question_order[block], "_", position)
    risk_sev_var <- paste0("risk_sev_", block, "_", question_order[block], "_", position)
    
    # Ensure these are numeric
    conjoint_rp[[risk_prob_var]] <- as.numeric(as.character(conjoint_rp[[risk_prob_var]]))
    conjoint_rp[[risk_sev_var]] <- as.numeric(as.character(conjoint_rp[[risk_sev_var]]))
  }
}

### Populate the new `rprob_*` and `rs_*` variables for flooding, infectious diseases, pesticides, and lifestyle
for (block in 1:24) {
  order <- strsplit(question_order[block], "")[[1]]  # Break the block into its components (f, i, p, l)
  
  for (position in 1:4) {
    risk_prob_var <- paste0("risk_prob_", block, "_", question_order[block], "_", position)
    risk_sev_var  <- paste0("risk_sev_", block, "_", question_order[block], "_", position)
    
    # Find indices where the value is not NA
    indices <- !is.na(conjoint_rp[[risk_prob_var]])
    
    # Assign values to the correct rprob_* and rs_* variables based on the type (f = flooding, i = infectious, etc.)
    if (order[position] == "f") {
      conjoint_rp$rprob_flooding[indices] <- conjoint_rp[[risk_prob_var]][indices]
      conjoint_rp$rs_flooding[indices] <- conjoint_rp[[risk_sev_var]][indices]
    } else if (order[position] == "i") {
      conjoint_rp$rprob_infectious[indices] <- conjoint_rp[[risk_prob_var]][indices]
      conjoint_rp$rs_infectious[indices] <- conjoint_rp[[risk_sev_var]][indices]
    } else if (order[position] == "p") {
      conjoint_rp$rprob_pesticides[indices] <- conjoint_rp[[risk_prob_var]][indices]
      conjoint_rp$rs_pesticides[indices] <- conjoint_rp[[risk_sev_var]][indices]
    } else if (order[position] == "l") {
      conjoint_rp$rprob_lifestyle[indices] <- conjoint_rp[[risk_prob_var]][indices]
      conjoint_rp$rs_lifestyle[indices] <- conjoint_rp[[risk_sev_var]][indices]
    }
  }
}


summary(conjoint_rp$rprob_flooding)
summary(conjoint_rp$rs_flooding)
summary(conjoint_rp$rprob_infectious)
summary(conjoint_rp$rs_infectious)
summary(conjoint_rp$rprob_pesticides)
summary(conjoint_rp$rs_pesticides)
summary(conjoint_rp$rprob_lifestyle)
summary(conjoint_rp$rs_lifestyle)


conjoint_rp <- conjoint_rp %>% drop_na(rprob_flooding)
conjoint_rp <- conjoint_rp %>% drop_na(rs_flooding)
conjoint_rp <- conjoint_rp %>% drop_na(rprob_infectious)
conjoint_rp <- conjoint_rp %>% drop_na(rs_infectious)
conjoint_rp <- conjoint_rp %>% drop_na(rprob_pesticides)
conjoint_rp <- conjoint_rp %>% drop_na(rs_pesticides)
conjoint_rp <- conjoint_rp %>% drop_na(rprob_lifestyle)
conjoint_rp <- conjoint_rp %>% drop_na(rs_lifestyle)


###Topic-specific risk perception variables
conjoint_rp$rp_flooding <- conjoint_rp$rprob_flooding * conjoint_rp$rs_flooding
conjoint_rp$rp_infectious <- conjoint_rp$rprob_infectious * conjoint_rp$rs_infectious
conjoint_rp$rp_pesticides <- conjoint_rp$rprob_pesticides * conjoint_rp$rs_pesticides
conjoint_rp$rp_lifestyle <- conjoint_rp$rprob_lifestyle * conjoint_rp$rs_lifestyle

summary(conjoint_rp$rp_flooding)
summary(conjoint_rp$rp_infectious)
summary(conjoint_rp$rp_pesticides)
summary(conjoint_rp$rp_lifestyle)


#RECODING OF RP VARIABLES
conjoint_rp <- conjoint_rp %>%
  mutate(
    rp_cat_flood = ifelse(rp_flooding >= 0 & rp_flooding < 6, "Low flood RP",
              ifelse(rp_flooding >= 6 & rp_flooding <= 16, "High flood RP", NA)),
    rp_cat_flood = factor(rp_cat_flood, levels = c("Low flood RP", "High flood RP"))  # Convert rp_cat to a factor
  )


conjoint_rp <- conjoint_rp %>%
  mutate(
    rp_cat_infect = ifelse(rp_infectious >= 0 & rp_infectious < 6, "Low infectious disease RP",
              ifelse(rp_infectious >= 6 & rp_infectious <= 16, "High infectious disease RP", NA)),
    rp_cat_infect = factor(rp_cat_infect, levels = c("Low infectious disease RP", "High infectious disease RP"))  # Convert rp_cat to a factor
  )

conjoint_rp <- conjoint_rp %>%
  mutate(
    rp_cat_pest = ifelse(rp_pesticides >= 0 & rp_pesticides < 6, "Low pesticide RP",
              ifelse(rp_pesticides >= 6 & rp_pesticides <= 16, "High pesticide RP", NA)),
    rp_cat_pest = factor(rp_cat_pest, levels = c("Low pesticide RP", "High pesticide RP"))  # Convert rp_cat to a factor
  )

conjoint_rp <- conjoint_rp %>%
  mutate(
    rp_cat_life = ifelse(rp_lifestyle >= 0 & rp_lifestyle < 6, "Low lifestyle disease RP",
              ifelse(rp_lifestyle >= 6 & rp_lifestyle <= 16, "High lifestyle disease RP", NA)),
    rp_cat_life = factor(rp_cat_life, levels = c("Low lifestyle disease RP", "High lifestyle disease RP"))  # Convert rp_cat to a factor
  )

summary(conjoint_rp$rp_cat_flood)
summary(conjoint_rp$rp_cat_infect)
summary(conjoint_rp$rp_cat_pest)
summary(conjoint_rp$rp_cat_life)

# Group by respondent and get the dominant category for each respondent
respondent_rp_flood <- conjoint_rp %>%
  group_by(Response.ID) %>%
  summarize(
    dominant_rp_cat_flood = ifelse(sum(rp_cat_flood == "High flood RP", na.rm = TRUE) > sum(rp_cat_flood == "Low flood RP", na.rm = TRUE), 
                             "High flood RP", 
                             "Low flood RP")
  )

respondent_rp_pest <- conjoint_rp %>%
  group_by(Response.ID) %>%
  summarize(
    dominant_rp_cat_pest = ifelse(sum(rp_cat_pest == "High pesticide RP", na.rm = TRUE) > sum(rp_cat_pest == "Low pesticide RP", na.rm = TRUE), 
                             "High pesticide RP", 
                             "Low pesticide RP")
  )

respondent_rp_infect <- conjoint_rp %>%
  group_by(Response.ID) %>%
  summarize(
    dominant_rp_cat_infect = ifelse(sum(rp_cat_infect == "High infectious disease RP", na.rm = TRUE) > sum(rp_cat_infect == "Low infectious disease RP", na.rm = TRUE), 
                             "High infectious disease RP", 
                             "Low infectious disease RP")
  )

respondent_rp_life <- conjoint_rp %>%
  group_by(Response.ID) %>%
  summarize(
    dominant_rp_cat_life = ifelse(sum(rp_cat_life == "High lifestyle disease RP", na.rm = TRUE) > sum(rp_cat_life == "Low lifestyle disease RP", na.rm = TRUE), 
                             "High lifestyle disease RP", 
                             "Low lifestyle disease RP")
  )

# Check the distribution of high and low risk perception by respondents
table(respondent_rp_flood$dominant_rp_cat_flood)
table(respondent_rp_pest$dominant_rp_cat_pest)
table(respondent_rp_infect$dominant_rp_cat_infect)
table(respondent_rp_life$dominant_rp_cat_life)

###GENERAL RISK PERCEPTION VARIABLE

# Calculate the average of the four composite variables
conjoint_rp$risk_perception <- rowMeans(conjoint_rp[, c("rp_flooding", 
                                                   "rp_infectious", 
                                                   "rp_pesticides", 
                                                   "rp_lifestyle")], na.rm = TRUE)

mean(conjoint_rp$risk_perception, na.rm = TRUE)

summary(conjoint_rp$risk_perception)

#The possible scores are as follows:
#"0" (completely 'improbable' regardless of severity and vice versa)
#"1" ('improbable' AND 'not severe')
#"2" ('improbable' AND 'neither severe not not severe', 'neither probable nor improbable' AND 'not severe')
#"3" ('improbable' AND 'severe', 'probable' AND 'not severe')
#"4" ('neither probable nor improbable' AND 'neither severe not not severe', 'very probable' AND 'not severe', 'not probable' AND 'severe')
#"6" ('probable' AND 'neither severe not not severe', 'neither probable nor improbable' AND 'severe')
#"8" ('very probable' AND 'neither severe not not severe', 'neither probable nor improbable' AND very 'severe')
#"9" ('probable' AND 'severe')
#"12"('very probable' AND 'severe', 'probable' AND 'very severe')
#"16"('very probable' AND 'very severe')
#The middle group is fairly large and makes it difficult to separate this variable into two categories. As such, two models are run: one with two categories and one with three.

conjoint_rp <- conjoint_rp %>%
  mutate(
    rp_cat = ifelse(risk_perception >= 0 & risk_perception < 6, "Low RP",
              ifelse(risk_perception >= 6 & risk_perception <= 16, "High RP", NA)),
    rp_cat = factor(rp_cat, levels = c("Low RP", "High RP"))  # Convert age_cat to a factor
  )

summary(conjoint_rp$rp_cat)


# Assuming 'respondent_id' is the unique identifier for each respondent
# Group by respondent and get the dominant category for each respondent
respondent_risk_perception <- conjoint_rp %>%
  group_by(Response.ID) %>%
  summarize(
    dominant_rp_cat = ifelse(sum(rp_cat == "High RP", na.rm = TRUE) > sum(rp_cat == "Low RP", na.rm = TRUE), 
                             "High RP", 
                             "Low RP")
  )

# Check the distribution of high and low risk perception by respondents
table(respondent_risk_perception$dominant_rp_cat)



```

#RENAME VARIABLES

```{r}

# Renaming attributes and levels in English
conjoint_rp <-  conjoint_rp %>% dplyr::rename(Topic=Onderwerp, 
                                            Relation.to.you=Relatie.tot.u,
                                            Apparent.knowledge.of.the.person=Schijnbare.kennis.van.de.persoon,
                                            Apparent.concern.of.the.person=Schijnbare.bezorgdheid.van.de.persoon,
                                            Apparent.motivation.of.the.person=Schijnbare.motivatie.van.de.persoon,
                                            Closeness.with.you=Band.met.u,
                                            Gender.of.the.person=Geslacht.van.de.persoon
                                            ) %>%
                                     mutate(
                                            Topic = recode_factor(Topic,
                                                  "overstromingen" = "floods", 
                                                  "besmettelijke ziekten" = "infectious diseases",
                                                  "pesticiden" = "pesticides",
                                                  "levensstijlziekten" = "lifestyle diseases",),
                                            Gender.of.the.person = recode_factor(Gender.of.the.person,
                                                  "vrouw" = "female", 
                                                  "man" = "male"),
                                            Relation.to.you = recode_factor(Relation.to.you,
                                                  "vriend/vriendin (platonisch)" = "friend",
                                                  "collega/studiegenoot" = "colleague/fellow student",
                                                  "familielid" = "relative"),
                                            Closeness.with.you = recode_factor(Closeness.with.you,
                                                  "zeer hecht" = "very close", 
                                                  "redelijk hecht" = "fairly close",
                                                  "enigszins hecht" = "somewhat close",
                                                  "niet hecht" = "not close"),
                                            Apparent.knowledge.of.the.person = recode_factor(Apparent.knowledge.of.the.person,
                                                  "zeer goed geinformeerd" = "very well-informed", 
                                                  "redelijk goed geinformeerd" = "fairly well-informed",
                                                  "enigszins geinformeerd" = "somewhat well-informed",
                                                  "niet geinformeerd" = "not well-informed"),
                                            Apparent.concern.of.the.person = recode_factor(Apparent.concern.of.the.person,
                                                  "zeer bezorgd" = "very concerned", 
                                                  "redelijk bezorgd" = "fairly concerned",
                                                  "enigszins bezorgd" = "somewhat concerned",
                                                  "niet bezorgd" = "not concerned"),
                                            Apparent.motivation.of.the.person = recode_factor(Apparent.motivation.of.the.person,
                                                  "u te overtuigen om meer maatregelen te nemen" = "to convince you to take more precautions", 
                                                  "informatie over het risico met u te willen uitwisselen" = "to exchange information with you",
                                                  "zich minder ongerust voelen" = "to feel less anxious",
                                                  "u te overtuigen om minder maatregelen te nemen" = "to convince you to take fewer precautions"))


#Check renamings
table(conjoint_rp$Topic)
table(conjoint_rp$Gender.of.the.person)
table(conjoint_rp$Relation.to.you)
table(conjoint_rp$Closeness.with.you)
table(conjoint_rp$Apparent.knowledge.of.the.person)
table(conjoint_rp$Apparent.concern.of.the.person)
table(conjoint_rp$Apparent.motivation.of.the.person)


```

#CHANGE VARIABLE ORDER AND SET REFERENCE CATEGORIES

```{r}
#REORDER LEVELS
conjoint_rp$Topic <-                             factor(conjoint_rp$Topic,
                                                  levels = c("pesticides", "floods", "lifestyle diseases", "infectious diseases"))
conjoint_rp$Gender.of.the.person <-              factor(conjoint_rp$Gender.of.the.person,
                                                  levels = c(  "male", "female"))
conjoint_rp$Relation.to.you <-                   factor(conjoint_rp$Relation.to.you,
                                                  levels = c("colleague/fellow student", "friend", "relative"))
conjoint_rp$Closeness.with.you <-                factor(conjoint_rp$Closeness.with.you,
                                                  levels = c("not close", "somewhat close", "fairly close", "very close"))
conjoint_rp$Apparent.knowledge.of.the.person <-  factor(conjoint_rp$Apparent.knowledge.of.the.person,
                                                  levels = c("not well-informed", "somewhat well-informed", "fairly well-informed", "very well-informed"))
conjoint_rp$Apparent.concern.of.the.person <-    factor(conjoint_rp$Apparent.concern.of.the.person,
                                                  levels = c("not concerned", "somewhat concerned", "fairly concerned", "very concerned"))
conjoint_rp$Apparent.motivation.of.the.person <- factor(conjoint_rp$Apparent.motivation.of.the.person,
                                                  levels = c("to convince you to take more precautions", "to convince you to take fewer precautions", "to feel less anxious", "to exchange information with you"))


```

######RUN MAIN AMCE, MM, AND DIFFERENCE IN MM (AVERAGE RP)

```{r}
###AMCE

#Full model
amce_rp <- cj(conjoint_rp, selected ~ 
             Topic +
             Gender.of.the.person +
             Relation.to.you +
             Closeness.with.you + 
             Apparent.knowledge.of.the.person +
             Apparent.concern.of.the.person +
             Apparent.motivation.of.the.person, 
                  by = ~ rp_cat, 
                  id = ~ Response.ID,
                  estimate = "amce")


# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

amce_rp_plot <- ggplot(amce_rp, aes(x = estimate, y = level, color = rp_cat, shape = rp_cat)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("darkolivegreen2", "darkgreen")) + 
  scale_shape_manual(values = c(16, 17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "AMCE", y = "All attributes by general risk perception", color = "General risk perception", shape = "General risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  ) +
  xlim(-0.3, 0.3)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/amce_rp.png", plot = amce_rp_plot, width = 8, height = 6, dpi = 300)


###MARGINAL MEANS

mm_rp <- cj(conjoint_rp, selected ~ 
                  Topic +
                  Gender.of.the.person +
                  Relation.to.you +
                  Closeness.with.you + 
                  Apparent.knowledge.of.the.person +
                  Apparent.concern.of.the.person +
                  Apparent.motivation.of.the.person, 
                  by = ~ rp_cat, 
                  id = ~ Response.ID,
                  estimate = "mm")


mm_rp$feature <- factor(mm_rp$feature,
                              levels = c("Topic",
                                         "Relation.to.you",
                                         "Apparent.knowledge.of.the.person",
                                         "Apparent.concern.of.the.person",
                                         "Apparent.motivation.of.the.person",
                                         "Closeness.with.you",
                                         "Gender.of.the.person"))

# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

mm_rp_plot <- ggplot(mm_rp, aes(x = estimate, y = level, color = rp_cat, shape = rp_cat)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("darkolivegreen2", "darkgreen")) +  # Custom colors for gender_cat
  scale_shape_manual(values = c(16, 17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "Marginal Means", y = "All attributes by general risk perception", color = "General risk perception", shape = "General risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  ) +
  xlim(0.3, 0.7)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/mm_rp.png", plot = mm_rp_plot, width = 8, height = 6, dpi = 300)


###DIFFERENCE IN MARGINAL MEANS

diff_rp <- cj(conjoint_rp, selected ~ 
                  Topic +
                  Gender.of.the.person +
                  Relation.to.you +
                  Closeness.with.you + 
                  Apparent.knowledge.of.the.person +
                  Apparent.concern.of.the.person +
                  Apparent.motivation.of.the.person, 
                  by = ~ rp_cat, 
                  id = ~ Response.ID,
                  estimate = "mm_differences")


diff_rp$feature <- factor(diff_rp$feature,
                              levels = c("Topic",
                                         "Relation.to.you",
                                         "Apparent.knowledge.of.the.person",
                                         "Apparent.concern.of.the.person",
                                         "Apparent.motivation.of.the.person",
                                         "Closeness.with.you",
                                         "Gender.of.the.person"))

# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

diff_rp_plot <- ggplot(diff_rp, aes(x = estimate, y = level, color = rp_cat, shape = rp_cat)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("darkgreen")) +  # Custom colors for gender_cat
  scale_shape_manual(values = c(17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "Differences in Marginal Means", y = "All attributes by general risk perception", color = "General risk perception", shape = "General risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  ) +
  xlim(-0.15, 0.15)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/diff_rp.png", plot = diff_rp_plot, width = 8, height = 6, dpi = 300)

```

######RUN MAIN AMCE, MM, AND DIFFERENCE IN MM (FLOOD RP)

```{r}
###AMCE

#Full model
amce_rp_flood <- cj(conjoint_rp, selected ~ 
             Topic +
             Gender.of.the.person +
             Relation.to.you +
             Closeness.with.you + 
             Apparent.knowledge.of.the.person +
             Apparent.concern.of.the.person +
             Apparent.motivation.of.the.person, 
                  by = ~ rp_cat_flood, 
                  id = ~ Response.ID,
                  estimate = "amce")


# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

amce_rp_flood_plot <- ggplot(amce_rp_flood, aes(x = estimate, y = level, color = rp_cat_flood, shape = rp_cat_flood)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("lightgreen", "darkolivegreen")) + 
  scale_shape_manual(values = c(16, 17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "AMCE", y = "All attributes by flood risk perception", color = "Flood risk perception", shape = "Flood risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  ) +
  xlim(-0.3, 0.3)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/amce_rp_flood.png", plot = amce_rp_flood_plot, width = 8, height = 6, dpi = 300)


###MARGINAL MEANS

mm_rp_flood <- cj(conjoint_rp, selected ~ 
                  Topic +
                  Gender.of.the.person +
                  Relation.to.you +
                  Closeness.with.you + 
                  Apparent.knowledge.of.the.person +
                  Apparent.concern.of.the.person +
                  Apparent.motivation.of.the.person, 
                  by = ~ rp_cat_flood, 
                  id = ~ Response.ID,
                  estimate = "mm")


mm_rp_flood$feature <- factor(mm_rp_flood$feature,
                              levels = c("Topic",
                                         "Relation.to.you",
                                         "Apparent.knowledge.of.the.person",
                                         "Apparent.concern.of.the.person",
                                         "Apparent.motivation.of.the.person",
                                         "Closeness.with.you",
                                         "Gender.of.the.person"))

# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

mm_rp_flood_plot <- ggplot(mm_rp_flood, aes(x = estimate, y = level, color = rp_cat_flood, shape = rp_cat_flood)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("lightgreen", "darkolivegreen")) + 
  scale_shape_manual(values = c(16, 17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "Marginal Means", y = "All attributes by flood risk perception", color = "Flood risk perception", shape = "Flood risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  ) +
  xlim(0.3, 0.7)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/mm_rp_flood.png", plot = mm_rp_flood_plot, width = 8, height = 6, dpi = 300)

###MARGINAL MEANS

diff_rp_flood <- cj(conjoint_rp, selected ~ 
                  Topic +
                  Gender.of.the.person +
                  Relation.to.you +
                  Closeness.with.you + 
                  Apparent.knowledge.of.the.person +
                  Apparent.concern.of.the.person +
                  Apparent.motivation.of.the.person, 
                  by = ~ rp_cat_flood, 
                  id = ~ Response.ID,
                  estimate = "mm_differences")


diff_rp_flood$feature <- factor(diff_rp_flood$feature,
                              levels = c("Topic",
                                         "Relation.to.you",
                                         "Apparent.knowledge.of.the.person",
                                         "Apparent.concern.of.the.person",
                                         "Apparent.motivation.of.the.person",
                                         "Closeness.with.you",
                                         "Gender.of.the.person"))

# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

diff_rp_flood_plot <- ggplot(diff_rp_flood, aes(x = estimate, y = level, color = rp_cat_flood, shape = rp_cat_flood)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("darkolivegreen")) + 
  scale_shape_manual(values = c(17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "Difference in Marginal Means", y = "All attributes by flood risk perception", color = "Flood risk perception", shape = "Flood risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  ) +
  xlim(-0.15, 0.15)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/diff_rp_flood.png", plot = diff_rp_flood_plot, width = 8, height = 6, dpi = 300)



```

######RUN MAIN AMCE, MM, AND DIFFERENCE IN MM (PESTICIDE RP)

```{r}
###AMCE

#Full model
amce_rp_pest <- cj(conjoint_rp, selected ~ 
             Topic +
             Gender.of.the.person +
             Relation.to.you +
             Closeness.with.you + 
             Apparent.knowledge.of.the.person +
             Apparent.concern.of.the.person +
             Apparent.motivation.of.the.person, 
                  by = ~ rp_cat_pest, 
                  id = ~ Response.ID,
                  estimate = "amce")


# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

amce_rp_pest_plot <- ggplot(amce_rp_pest, aes(x = estimate, y = level, color = rp_cat_pest, shape = rp_cat_pest)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("lightgreen", "darkolivegreen")) + 
  scale_shape_manual(values = c(16, 17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "AMCE", y = "All attributes by pesticide risk perception", color = "Pesticide risk perception", shape = "Pesticide risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  ) +
  xlim(-0.3, 0.3)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/amce_rp_pest.png", plot = amce_rp_pest_plot, width = 8, height = 6, dpi = 300)


###MARGINAL MEANS

mm_rp_pest <- cj(conjoint_rp, selected ~ 
                  Topic +
                  Gender.of.the.person +
                  Relation.to.you +
                  Closeness.with.you + 
                  Apparent.knowledge.of.the.person +
                  Apparent.concern.of.the.person +
                  Apparent.motivation.of.the.person, 
                  by = ~ rp_cat_pest, 
                  id = ~ Response.ID,
                  estimate = "mm")


mm_rp_pest$feature <- factor(mm_rp_pest$feature,
                              levels = c("Topic",
                                         "Relation.to.you",
                                         "Apparent.knowledge.of.the.person",
                                         "Apparent.concern.of.the.person",
                                         "Apparent.motivation.of.the.person",
                                         "Closeness.with.you",
                                         "Gender.of.the.person"))

# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

mm_rp_pest_plot <- ggplot(mm_rp_pest, aes(x = estimate, y = level, color = rp_cat_pest, shape = rp_cat_pest)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("lightgreen", "darkolivegreen")) + 
  scale_shape_manual(values = c(16, 17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "Marginal Means", y = "All attributes by pesticide risk perception", color = "Pesticide risk perception", shape = "Pesticide risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  ) +
  xlim(0.3, 0.7)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/mm_rp_pest.png", plot = mm_rp_pest_plot, width = 8, height = 6, dpi = 300)


###DIFFERENCE IN MARGINAL MEANS

diff_rp_pest <- cj(conjoint_rp, selected ~ 
                  Topic +
                  Gender.of.the.person +
                  Relation.to.you +
                  Closeness.with.you + 
                  Apparent.knowledge.of.the.person +
                  Apparent.concern.of.the.person +
                  Apparent.motivation.of.the.person, 
                  by = ~ rp_cat_pest, 
                  id = ~ Response.ID,
                  estimate = "mm_differences")


diff_rp_pest$feature <- factor(diff_rp_pest$feature,
                              levels = c("Topic",
                                         "Relation.to.you",
                                         "Apparent.knowledge.of.the.person",
                                         "Apparent.concern.of.the.person",
                                         "Apparent.motivation.of.the.person",
                                         "Closeness.with.you",
                                         "Gender.of.the.person"))

# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

diff_rp_pest_plot <- ggplot(diff_rp_pest, aes(x = estimate, y = level, color = rp_cat_pest, shape = rp_cat_pest)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("darkolivegreen")) + 
  scale_shape_manual(values = c(17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "Difference in Marginal Means", y = "All attributes by pesticide risk perception", color = "Pesticide risk perception", shape = "Pesticide risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  ) +
  xlim(-0.15, 0.15)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/diff_rp_pest.png", plot = diff_rp_pest_plot, width = 8, height = 6, dpi = 300)

```

######RUN MAIN AMCE AND MM (INFECTIOUS DISEASE RP)

```{r}
###AMCE

#Full model
amce_rp_infect <- cj(conjoint_rp, selected ~ 
             Topic +
             Gender.of.the.person +
             Relation.to.you +
             Closeness.with.you + 
             Apparent.knowledge.of.the.person +
             Apparent.concern.of.the.person +
             Apparent.motivation.of.the.person, 
                  by = ~ rp_cat_infect, 
                  id = ~ Response.ID,
                  estimate = "amce")


# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

amce_rp_infect_plot <- ggplot(amce_rp_infect, aes(x = estimate, y = level, color = rp_cat_infect, shape = rp_cat_infect)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("lightgreen", "darkolivegreen")) + 
  scale_shape_manual(values = c(16, 17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "AMCE", y = "All attributes by infectious disease risk perception", color = "Infectious disease risk perception", shape = "Infectious disease risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  ) +
  xlim(-0.3, 0.3)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/amce_rp_infect.png", plot = amce_rp_infect_plot, width = 8, height = 6, dpi = 300)


###MARGINAL MEANS

mm_rp_infect <- cj(conjoint_rp, selected ~ 
                  Topic +
                  Gender.of.the.person +
                  Relation.to.you +
                  Closeness.with.you + 
                  Apparent.knowledge.of.the.person +
                  Apparent.concern.of.the.person +
                  Apparent.motivation.of.the.person, 
                  by = ~ rp_cat_infect, 
                  id = ~ Response.ID,
                  estimate = "mm")


mm_rp_infect$feature <- factor(mm_rp_infect$feature,
                              levels = c("Topic",
                                         "Relation.to.you",
                                         "Apparent.knowledge.of.the.person",
                                         "Apparent.concern.of.the.person",
                                         "Apparent.motivation.of.the.person",
                                         "Closeness.with.you",
                                         "Gender.of.the.person"))

# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

mm_rp_infect_plot <- ggplot(mm_rp_infect, aes(x = estimate, y = level, color = rp_cat_infect, shape = rp_cat_infect)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("lightgreen", "darkolivegreen")) + 
  scale_shape_manual(values = c(16, 17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "Marginal Means", y = "All attributes by infectious disease risk perception", color = "Infectious disease risk perception", shape = "Infectious disease risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  )+
  xlim(0.3, 0.7)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/mm_rp_infect.png", plot = mm_rp_infect_plot, width = 8, height = 6, dpi = 300)

###MARGINAL MEANS

diff_rp_infect <- cj(conjoint_rp, selected ~ 
                  Topic +
                  Gender.of.the.person +
                  Relation.to.you +
                  Closeness.with.you + 
                  Apparent.knowledge.of.the.person +
                  Apparent.concern.of.the.person +
                  Apparent.motivation.of.the.person, 
                  by = ~ rp_cat_infect, 
                  id = ~ Response.ID,
                  estimate = "mm_difference")


diff_rp_infect$feature <- factor(diff_rp_infect$feature,
                              levels = c("Topic",
                                         "Relation.to.you",
                                         "Apparent.knowledge.of.the.person",
                                         "Apparent.concern.of.the.person",
                                         "Apparent.motivation.of.the.person",
                                         "Closeness.with.you",
                                         "Gender.of.the.person"))

# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

diff_rp_infect_plot <- ggplot(diff_rp_infect, aes(x = estimate, y = level, color = rp_cat_infect, shape = rp_cat_infect)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("darkolivegreen")) + 
  scale_shape_manual(values = c(17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "Difference in Marginal Means", y = "All attributes by infectious disease risk perception", color = "Infectious disease risk perception", shape = "Infectious disease risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  )+
  xlim(-0.15, 0.15)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/diff_rp_infect.png", plot = diff_rp_infect_plot, width = 8, height = 6, dpi = 300)

```

######RUN MAIN AMCE, MM, AND DIFFERENCE IN MM (LIFESTYLE DISEASE RP)

```{r}
###AMCE

#Full model
amce_rp_life <- cj(conjoint_rp, selected ~ 
             Topic +
             Gender.of.the.person +
             Relation.to.you +
             Closeness.with.you + 
             Apparent.knowledge.of.the.person +
             Apparent.concern.of.the.person +
             Apparent.motivation.of.the.person, 
                  by = ~ rp_cat_life, 
                  id = ~ Response.ID,
                  estimate = "amce")


# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

amce_rp_life_plot <- ggplot(amce_rp_life, aes(x = estimate, y = level, color = rp_cat_life, shape = rp_cat_life)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("lightgreen", "darkolivegreen")) + 
  scale_shape_manual(values = c(16, 17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "AMCE", y = "All attributes by lifestyle disease risk perception", color = "Lifestyle disease risk perception", shape = "Lifestyle disease risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  ) +
  xlim(-0.3, 0.3)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/amce_rp_life.png", plot = amce_rp_life_plot, width = 8, height = 6, dpi = 300)


###MARGINAL MEANS

mm_rp_life <- cj(conjoint_rp, selected ~ 
                  Topic +
                  Gender.of.the.person +
                  Relation.to.you +
                  Closeness.with.you + 
                  Apparent.knowledge.of.the.person +
                  Apparent.concern.of.the.person +
                  Apparent.motivation.of.the.person, 
                  by = ~ rp_cat_life, 
                  id = ~ Response.ID,
                  estimate = "mm")


mm_rp_life$feature <- factor(mm_rp_life$feature,
                              levels = c("Topic",
                                         "Relation.to.you",
                                         "Apparent.knowledge.of.the.person",
                                         "Apparent.concern.of.the.person",
                                         "Apparent.motivation.of.the.person",
                                         "Closeness.with.you",
                                         "Gender.of.the.person"))

# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

mm_rp_life_plot <- ggplot(mm_rp_life, aes(x = estimate, y = level, color = rp_cat_life, shape = rp_cat_life)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c("lightgreen", "darkolivegreen")) + 
  scale_shape_manual(values = c(16, 17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "Marginal Means", y = "All attributes by lifestyle disease risk perception", color = "Lifestyle disease risk perception", shape = "Lifestyle disease risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  ) +
  xlim(0.3, 0.7)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/mm_rp_life.png", plot = mm_rp_life_plot, width = 8, height = 6, dpi = 300)


###DIFFERENCE IN MARGINAL MEANS

diff_rp_life <- cj(conjoint_rp, selected ~ 
                  Topic +
                  Gender.of.the.person +
                  Relation.to.you +
                  Closeness.with.you + 
                  Apparent.knowledge.of.the.person +
                  Apparent.concern.of.the.person +
                  Apparent.motivation.of.the.person, 
                  by = ~ rp_cat_life, 
                  id = ~ Response.ID,
                  estimate = "mm_difference")


diff_rp_life$feature <- factor(diff_rp_life$feature,
                              levels = c("Topic",
                                         "Relation.to.you",
                                         "Apparent.knowledge.of.the.person",
                                         "Apparent.concern.of.the.person",
                                         "Apparent.motivation.of.the.person",
                                         "Closeness.with.you",
                                         "Gender.of.the.person"))

# Define the dodge width to create separation between low and high risk perception lines
dodge_width <- 1

diff_rp_life_plot <- ggplot(diff_rp_life, aes(x = estimate, y = level, color = rp_cat_life, shape = rp_cat_life)) +
  geom_point(size = 2, position = position_dodge(width = dodge_width)) +  # Offset points
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 1, 
                 position = position_dodge(width = dodge_width)) +  # Offset error bars
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +  # Add a dashed vertical line at x = 0
  scale_color_manual(values = c( "darkolivegreen")) + 
  scale_shape_manual(values = c(17)) +  # Custom shapes: 16 = filled circle, 17 = filled triangle
  theme_minimal() +  # Apply a minimal theme
  labs(x = "Difference in Marginal Means", y = "All attributes by lifestyle disease risk perception", color = "Lifestyle disease risk perception", shape = "Lifestyle disease risk perception") + # Customize labels
  theme(
    panel.background = element_rect(fill = "white", color = NA),  # Set the background of the plot panel to white
    plot.background = element_rect(fill = "white", color = NA)  # Set the overall background of the plot to white
  ) +
  xlim(-0.15, 0.15)  # Set the x-axis limits (replace with your desired values)

ggsave("figures/risk_perception/diff_rp_life.png", plot = diff_rp_life_plot, width = 8, height = 6, dpi = 300)

```

#BIVARIATE CORRELATIONS OF TOPIC-SPECIFIC RISK PERCEPTIONS

```{r}

# Create a dataframe with the five variables
risk_vars <- conjoint_rp[, c("risk_perception", "rp_flooding", "rp_infectious", "rp_pesticides", "rp_lifestyle")]

# Calculate the correlation matrix
correlation_matrix <- cor(risk_vars, use = "complete.obs", method = "pearson")

# Display the correlation matrix
print(correlation_matrix)


```

#####DIAGNOSTICS

```{r}

# BALANCE TESTING: GENERAL

# Manually map rp_cat to numeric
conjoint_rp$rp_cat_num <- ifelse(conjoint_rp$rp_cat == "Low RP", 0, 
                                      ifelse(conjoint_rp$rp_cat == "High RP", 1, NA))

# Check for any NA values introduced by this process
table(conjoint_rp$rp_cat_num)

# Calculate marginal means
rpbalance <- plot(mm(conjoint_rp, rp_cat_num ~ 
                 Topic +
                 Gender.of.the.person +
                 Apparent.motivation.of.the.person +
                 Apparent.concern.of.the.person +
                 Relation.to.you +
                 Closeness.with.you + 
                 Apparent.knowledge.of.the.person,
                 id = ~Response.ID), 
                xlim = c(0, 0.5), vline = mean(conjoint_rp$rp_cat_num, na.rm = TRUE))

# Adjust the plot for larger text and dots
rpbalance_adj <- rpbalance +
  geom_point(size = 3, position = position_dodge(width = 0.5)) +
    ggtitle("Balance Testing: General Risk Perception") +  # Add the title here
    labs(color = "Attribute") +  # Change the label for the color legend
  theme(
    text = element_text(size = 20),            # Increase base text size
    axis.title = element_text(size = 22),      # Increase axis title size
    axis.text = element_text(size = 18),       # Increase axis tick labels size
    legend.title = element_text(size = 22),    # Increase legend title size
    legend.text = element_text(size = 18),     # Increase legend text size
    plot.title = element_text(size = 24, face = "bold"),  # Increase plot title size and style
    strip.text = element_text(size = 18)       # Increase facet labels text size (if using facets)
  ) +
  guides(color = guide_legend(nrow = 4))  # Adjust the number of rows in the legend

# Save the plot
ggsave("diagnostics/balance_rp.png", plot = rpbalance_adj, width = 15, height = 11, dpi = 300)


# BALANCE TESTING: FLOOD

# Manually map rp_cat_flood to numeric
conjoint_rp$rp_cat_flood_num <- ifelse(conjoint_rp$rp_cat_flood == "Low flood RP", 0, 
                                      ifelse(conjoint_rp$rp_cat_flood == "High flood RP", 1, NA))

# Check for any NA values introduced by this process
table(conjoint_rp$rp_cat_flood_num)

# Calculate marginal means
rpbalance_flood <- plot(mm(conjoint_rp, rp_cat_flood_num ~ 
                 Topic +
                 Gender.of.the.person +
                 Apparent.motivation.of.the.person +
                 Apparent.concern.of.the.person +
                 Relation.to.you +
                 Closeness.with.you + 
                 Apparent.knowledge.of.the.person,
                 id = ~Response.ID), 
                xlim = c(0, 0.7), vline = mean(conjoint_rp$rp_cat_flood_num, na.rm = TRUE))

# Adjust the plot for larger text and dots
rpbalance_flood_adj <- rpbalance_flood +
  geom_point(size = 3, position = position_dodge(width = 0.5)) +
    ggtitle("Balance Testing: Flood Risk Perception") +  # Add the title here
    labs(color = "Attribute") +  # Change the label for the color legend
  theme(
    text = element_text(size = 20),            # Increase base text size
    axis.title = element_text(size = 22),      # Increase axis title size
    axis.text = element_text(size = 18),       # Increase axis tick labels size
    legend.title = element_text(size = 22),    # Increase legend title size
    legend.text = element_text(size = 18),     # Increase legend text size
    plot.title = element_text(size = 24, face = "bold"),  # Increase plot title size and style
    strip.text = element_text(size = 18)       # Increase facet labels text size (if using facets)
  ) +
  guides(color = guide_legend(nrow = 4))  # Adjust the number of rows in the legend

# Save the plot
ggsave("diagnostics/balance_rp_flood.png", plot = rpbalance_flood_adj, width = 15, height = 11, dpi = 300)



# BALANCE TESTING: PESTICIDES

# Manually map rp_cat_pest to numeric
conjoint_rp$rp_cat_pest_num <- ifelse(conjoint_rp$rp_cat_pest == "Low pesticide RP", 0, 
                                      ifelse(conjoint_rp$rp_cat_pest == "High pesticide RP", 1, NA))

# Check for any NA values introduced by this process
table(conjoint_rp$rp_cat_pest_num)

# Calculate marginal means
rpbalance_pest <- plot(mm(conjoint_rp, rp_cat_pest_num ~ 
                 Topic +
                 Gender.of.the.person +
                 Apparent.motivation.of.the.person +
                 Apparent.concern.of.the.person +
                 Relation.to.you +
                 Closeness.with.you + 
                 Apparent.knowledge.of.the.person,
                 id = ~Response.ID), 
                xlim = c(0, 0.7), vline = mean(conjoint_rp$rp_cat_pest_num, na.rm = TRUE))

# Adjust the plot for larger text and dots
rpbalance_pest_adj <- rpbalance_pest +
  geom_point(size = 3, position = position_dodge(width = 0.5)) +
    ggtitle("Balance Testing: Pesticide Risk Perception") +  # Add the title here
    labs(color = "Attribute") +  # Change the label for the color legend
  theme(
    text = element_text(size = 20),            # Increase base text size
    axis.title = element_text(size = 22),      # Increase axis title size
    axis.text = element_text(size = 18),       # Increase axis tick labels size
    legend.title = element_text(size = 22),    # Increase legend title size
    legend.text = element_text(size = 18),     # Increase legend text size
    plot.title = element_text(size = 24, face = "bold"),  # Increase plot title size and style
    strip.text = element_text(size = 18)       # Increase facet labels text size (if using facets)
  ) +
  guides(color = guide_legend(nrow = 4))  # Adjust the number of rows in the legend

# Save the plot
ggsave("diagnostics/balance_rp_pest.png", plot = rpbalance_pest_adj, width = 15, height = 11, dpi = 300)


# BALANCE TESTING: INFECTIOUS DISEASES

# Manually map rp_cat_infect to numeric
conjoint_rp$rp_cat_infect_num <- ifelse(conjoint_rp$rp_cat_infect == "Low infectious disease RP", 0, 
                                      ifelse(conjoint_rp$rp_cat_infect == "High infectious disease RP", 1, NA))

# Check for any NA values introduced by this process
table(conjoint_rp$rp_cat_infect_num)

# Calculate marginal means
rpbalance_infect <- plot(mm(conjoint_rp, rp_cat_infect_num ~ 
                 Topic +
                 Gender.of.the.person +
                 Apparent.motivation.of.the.person +
                 Apparent.concern.of.the.person +
                 Relation.to.you +
                 Closeness.with.you + 
                 Apparent.knowledge.of.the.person,
                 id = ~Response.ID), 
                xlim = c(0, 1), vline = mean(conjoint_rp$rp_cat_infect_num, na.rm = TRUE))

# Adjust the plot for larger text and dots
rpbalance_infect_adj <- rpbalance_infect +
  geom_point(size = 3, position = position_dodge(width = 0.5)) +
    ggtitle("Balance Testing: Infectious Disease Risk Perception") +  # Add the title here
    labs(color = "Attribute") +  # Change the label for the color legend
  theme(
    text = element_text(size = 20),            # Increase base text size
    axis.title = element_text(size = 22),      # Increase axis title size
    axis.text = element_text(size = 18),       # Increase axis tick labels size
    legend.title = element_text(size = 22),    # Increase legend title size
    legend.text = element_text(size = 18),     # Increase legend text size
    plot.title = element_text(size = 24, face = "bold"),  # Increase plot title size and style
    strip.text = element_text(size = 18)       # Increase facet labels text size (if using facets)
  ) +
  guides(color = guide_legend(nrow = 4))  # Adjust the number of rows in the legend

# Save the plot
ggsave("diagnostics/balance_rp_infect.png", plot = rpbalance_infect_adj, width = 15, height = 11, dpi = 300)


# BALANCE TESTING: LIFESTYLE DISEASES

# Manually map rp_cat_life to numeric
conjoint_rp$rp_cat_life_num <- ifelse(conjoint_rp$rp_cat_life == "Low lifestyle disease RP", 0, 
                                      ifelse(conjoint_rp$rp_cat_life == "High lifestyle disease RP", 1, NA))

# Check for any NA values introduced by this process
table(conjoint_rp$rp_cat_life_num)

# Calculate marginal means
rpbalance_life <- plot(mm(conjoint_rp, rp_cat_life_num ~ 
                 Topic +
                 Gender.of.the.person +
                 Apparent.motivation.of.the.person +
                 Apparent.concern.of.the.person +
                 Relation.to.you +
                 Closeness.with.you + 
                 Apparent.knowledge.of.the.person,
                 id = ~Response.ID), 
                xlim = c(0, 0.7), vline = mean(conjoint_rp$rp_cat_life_num, na.rm = TRUE))

# Adjust the plot for larger text and dots
rpbalance_life_adj <- rpbalance_life +
  geom_point(size = 3, position = position_dodge(width = 0.5)) +
    ggtitle("Balance Testing: Lifestyle Disease Risk Perception") +  # Add the title here
    labs(color = "Attribute") +  # Change the label for the color legend
  theme(
    text = element_text(size = 20),            # Increase base text size
    axis.title = element_text(size = 22),      # Increase axis title size
    axis.text = element_text(size = 18),       # Increase axis tick labels size
    legend.title = element_text(size = 22),    # Increase legend title size
    legend.text = element_text(size = 18),     # Increase legend text size
    plot.title = element_text(size = 24, face = "bold"),  # Increase plot title size and style
    strip.text = element_text(size = 18)       # Increase facet labels text size (if using facets)
  ) +
  guides(color = guide_legend(nrow = 4))  # Adjust the number of rows in the legend

# Save the plot
ggsave("diagnostics/balance_rp_life.png", plot = rpbalance_life_adj, width = 15, height = 11, dpi = 300)

```


